/*
 * __Simbox.h
 *
 *  Created on:2015-1-27
 *      Author: movin
 *  This file will define the simulation box class
 */
// This class provide a date with year, month and day
// you can use = to assign for it.
#ifndef ____Simbox
#define ____Simbox
#include <iostream>
#include <string>
#include "fields.h"
//#include "particles.h"
#include "cellinf.h"
#include "common.h"
#include "type.h"
#include "userdef.h"
///Class data used to store simulation time, to remember you when you had simulated, this is comtemprarily unavailable
class __Date{
    public:
	int year;
	int month;
	int day;
	__Date(int y = 2015, int m = 1, int d = 1){ year = y; month = m; day = d;}
	__Date operator=(const __Date& other){
		this->year = other.year;
		this->month = other.month;
		this->day = other.day;
		return *this;
	}
	~__Date(){}
};


ostream& operator<<(ostream& out, const __Date& mydate){
///输出重载，日期可以直接使用cout输出了
	out<<mydate.year<<"-"<<mydate.month<<"-"<<mydate.day;
	return out;
}

istream& operator>>(istream& in, __Date& mydate){
///输入重载，可以通过输入三个整数来给一个日期赋值
	in>>mydate.year>>mydate.month>>mydate.day;
	return in;
}


///初始化一个模拟盒子，这个模拟盒子包含cell信息，粒子信息，场信息。粒子种类信息也在这里。
class __Simbox{
public:
    ///这个名字将会输出数据文件名的头部
 	string simname;
    ///日期，目前没用
	__Date simdate;
    ///所有总模拟盒子的边界
    __Vect2<double> basexymax;
    ///所有总模拟盒子的边界
    __Vect2<double> basexymin;
    ///盒子的网格标号最大和最小
	__Vect2<int> xygridmax;
    ///盒子的网格标号最大和最小
	__Vect2<int> xygridmin;
    ///盒子中网格的数目
	__Vect2<int> xycell;
    ///盒子x的范围
	__Vect2<double> xrange;
    ///盒子y的范围
	__Vect2<double> yrange;
    ///盒子，xy最大值
	__Vect2<double> xymax;
    ///盒子xy最小值
	__Vect2<double> xymin;
    ///模拟时间开始和结束点
    __Vect2<double> time;
    ///这个position本来是准备在2d MPI 中记录box位置的，在一维domain分解中没有必要
    __Vect2<int> position;
    ///网格长宽
    __Vect2<double> dxy;
    ///时间步数，不能自己指定，而是根据给定的网格宽度和时间长度、Corant 条件算出。
	int t_steps;
    ///粒子种类，4为默认，但可以不给粒子小于4时，最后面的粒子种类将不会启用 0：电子，1：离子，2：正电子，3：光子
	int species;
    /// speices struct info
    vector<Particles> spptcl;
    ///dt 根据网格宽度和Corant 条件算出
	double delta_t;
    ///代表当前box初始化时粒子id最大值
	int maxid;
    ///种类数组，当数组中对应位数为0时，该种类粒子不启用，为1时启用
    vector<int> ppc;
    ///当flag = 0时，模拟不会启动，因此必须设置该位
    bool flag; // 判断是否有粒子
    ///电磁场边界类型，为下、上、左、右的顺序，0：反射边界，1：吸收边界，2：周期
    int bd_type[4];
    ///damplen 开始的网格位置，同样为下、上、左、右
    int startindex[4];
    ///damplen 存储场的吸收边界情况下四个边界的damp 网格长度， 顺序为下、上、左、右
    int damplen[4];
    /// 单个网格单种粒子最大数目
    int max_ppc;
    /// 允许合并粒子的能量差异
    double delta_e;
    /// 允许合并粒子的速度差异 ratio
    double delta_v;
    /// 0 不需要开启，１需要开启, 顺序0-syn, 1-brem, 2-bwpair, 3-bhpair
    int monte_flag[4];
    /// 可能有2个表0,eta, 1, chi, 2 f, 3 h(eta), 4 T(chi), 5 chi for T, 6 eta for h
    vector<vector<double> > tables1d;
    /// 0 则开启电磁相互作用, 1 不开启, 默认开启
    int disable_em;
    /// 2d 目前只保存P(eta, kappa), P(f, kappa)
    vector<vector<vector<double> > > tables2d;
    /// 保存韧致辐射和ｂｈ过程的能量，以及能量比值
    vector<vector<double> > bh_tables1d;
    /// 保存韧致辐射, bh截面
    vector<vector<vector<double> > > bh_tables2d;
    ///存储越界粒子，当前box中所有越过母网格的粒子，是个vector，包含所有类别的粒子
    vector<vector<__Ptcls> >crossing_bd_ptcls;
    /// filenames1d, filenames2d
    vector<string> filenames1d;
    vector<string> filenames2d;
    ///当前box中的所有网格的二维数组，通过指标可以直接访问
    __Cellgroup cellgroup;
    ///当前box中所有grid上的场，是fields的一个二维数组
    __Fieldgroup fieldgroup;
    /// 高阶电磁场插值到粒子
    void weight_to_ptcls();
	__Simbox(string& simuname, __Date& today, \
			__Vect2<int> xygridmin, __Vect2<int> xygridmax,\
            __Vect2<int> xycell, __Vect2<double> xrange, \
			__Vect2<double> yrange, __Vect2<double> time,\
			int& tsteps, double& deltat, const int& species, const vector<Particles>& spptcl, int maxid, \
            const vector<int>& ppcref, __Vect2<int> positionref, int disable_em = 0):simname(simuname),\
                                                                                     simdate(today),\
			xygridmin(xygridmin), xygridmax(xygridmax), xycell(xycell), xrange(xrange),\
            yrange(yrange), time(time), t_steps(tsteps), delta_t(deltat), position(positionref)\
            , disable_em(disable_em), maxid(maxid), species(species)
            {
                xymin.member[0] = xrange.member[0];
                xymin.member[1] = yrange.member[0];
                xymax.member[0] = xrange.member[1];
                xymax.member[1] = yrange.member[1];
                dxy.member[0] = (xrange.member[1] - xrange.member[0]) / xycell.member[0];
                dxy.member[1] = (yrange.member[1] - yrange.member[0]) / xycell.member[1];
                this->spptcl = spptcl;
                ppc = ppcref;
            }
    __Simbox(){}
    ///利用构造函数提供的数据进行初始化，初始化cell，粒子，和场
    void init();
    /// 0 代表读取成功
    int read_tables1d();
    /// 0 代表读取成功
    int read_tables2d();
    ///通过这个函数得到某一左边，sp类粒子的密度，从而进行初始化
    double get_density(const __Vect2<double>& xycoords, int sp);
    ///对越界粒子进行处理
    void migration(const int& myrank);
    ///从坐标到网格的转换
    __Vect2<int> coords_to_cell(const __Vect2<double>& xycoords);
    ///场拷贝到网格
    void field_copy_to_cell();
    ///网格的流拷贝到场
    void cell_copy_to_field(int domain, const int& myrank);
    ///进行演化
    void advancest1(const int& rank);
    void advancest2(const int& rank);
    void advancest3(const int& rank);
    void advancest4(const int& rank, const int& n);
    void advancest5(const int& rank);
    ///数据dump，目前没用到
    void dump();
    ///得到当前时间点sp类粒子的总数目
    int get_sp_num(int sp);
    ///初始化flag，必须调用
    void init_flag();
    ///设置场的边界条件
    void set_field_boundary(const int bd_typer[], const int startindexr[], const int damplenr[]);
    ///拷贝越界粒子数据到MPI数据
    void copy_to_pline(__Pline& plinel, __Pline& pliner, const int direction);
    ///拷贝ghostcell场数据到MPI数据
    void copy_to_fline(__Fline& flinel, __Fline& fliner, const int direction);
    ///拷贝ghostcell数据到MPI数据
    void copy_to_jline(__Jline& jlinel, __Jline& jliner, const int direction);
    ///拷贝MPI数据到边界场
    void fline_copy_to(__Fline& flinel, __Fline& fliner, const int direction);
    ///拷贝MPI数据到当前越界粒子束中
    void pline_copy_to(__Pline& plinel, __Pline& pliner, const int direction);
    ///拷贝MPI网格数据到当前cellgroup数据中
    void jline_copy_to(__Jline& jlinel, __Jline& jliner, const int direction);
    /// 设置输入数据文件名
    void set_data_file(const vector<string>& fnames1d, const vector<string>& fnames2d);
    /// 拷贝ｆｌａｇ
    void set_monte_carlo(int monte_flags[])
    {
        memcpy(monte_flag, monte_flags, sizeof(int) * 4);
    }
    ///对内存进行优化
    void squeeze();
    /// 对sp类型的粒子进行合并
    void merge(const int& sp){cellgroup.merge_ptcls(max_ppc, delta_e, delta_v, sp); }
    void wide_merge(const int& sp, const int& cellnum);
	~__Simbox(){}
};


void __Simbox::set_data_file(const vector<string>& fnames1d, const vector<string>& fnames2d)
{
    for(int i = 0; i < fnames1d.size(); i ++)
    {
        filenames1d.push_back(fnames1d[i]);
    }
    for(int i = 0; i < fnames2d.size(); i ++)
    {
        filenames2d.push_back(fnames2d[i]);
    }
}


void __Simbox::set_field_boundary(const int bd_typer[], const int startindexr[], const int damplenr[])
{
    memcpy(bd_type, bd_typer, sizeof(int) * 4);
    memcpy(startindex, startindexr, sizeof(int) * 4);
    memcpy(damplen, damplenr, sizeof(int) * 4);
    memcpy(fieldgroup.bd_type, bd_typer, sizeof(int) * 4);
    memcpy(fieldgroup.startindex, startindexr, sizeof(int) * 4);
    memcpy(fieldgroup.damp_len, damplenr, sizeof(int) * 4);
}

void __Simbox::squeeze()
{
    cellgroup.squeeze();
    for(int i = 0; i < crossing_bd_ptcls.size(); i ++)
    {
        vector<__Ptcls>(crossing_bd_ptcls[i]).swap(crossing_bd_ptcls[i]);
    }
}

void __Simbox::init_flag()
{
    flag = (get_sp_num(0) || get_sp_num(1) || get_sp_num(2) || get_sp_num(3) != 0);
}

int __Simbox::get_sp_num(int sp)
{
    int retnum = 0;
    for(int i = 1; i < cellgroup.member.size() - 1; i ++)
    {
        for(int j = 1; j < cellgroup.member[i].size() - 1; j ++)
        {
            retnum += cellgroup.member[i][j].plist[sp].size();
            //cout<<" size "<<i<<"\t"<<j<<"\t"<<cellgroup.member[i][j].plist[sp].size()<<endl;
        }
    }

    return retnum;
}

double __Simbox::get_density(const __Vect2<double>& xycoords, int sp)
{
    double density  = 0;
    return density;
}

/// domain 代表是中心 0 还是边界 1
void __Simbox::cell_copy_to_field(int domain, const int& myrank)
{
    for(int i = 0; i < fieldgroup.member.size(); i ++)
    {
	    for(int j = 0; j < fieldgroup.member[i].size(); j ++)
	    {
            fieldgroup.member[i][j].jfield = cellgroup.jfield[i][j];
	        fieldgroup.member[i][j].rho = cellgroup.rhofield[i][j];
        }
    }
}



void __Simbox::field_copy_to_cell()
{
    for(int i = 0; i < cellgroup.member.size(); i ++){
        for(int j = 0; j < cellgroup.member[i].size(); j ++){
            cellgroup.member[i][j].fields.exfield.member[0] = fieldgroup.member[i][j].efield.\
                                                              member[0];
            cellgroup.member[i][j].fields.exfield.member[1] = fieldgroup.member[i][j+1].efield.\
                                                              member[0];
            cellgroup.member[i][j].fields.eyfield.member[0] = fieldgroup.member[i][j].efield.\
                                                              member[1];
            cellgroup.member[i][j].fields.eyfield.member[1] = fieldgroup.member[i+1][j].efield.\
                                                              member[1];
            cellgroup.member[i][j].fields.ezfield = fieldgroup.member[i][j].efield.member[2];
            cellgroup.member[i][j].fields.bxfield.member[0] = fieldgroup.member[i][j].bfield.\
                                                              member[0];
            cellgroup.member[i][j].fields.bxfield.member[1] = fieldgroup.member[i][j+1].bfield.\
                                                              member[0];
            cellgroup.member[i][j].fields.byfield.member[0] = fieldgroup.member[i][j].bfield.\
                                                              member[1];
            cellgroup.member[i][j].fields.byfield.member[1] = fieldgroup.member[i+1][j].bfield.\
                                                              member[1];
            cellgroup.member[i][j].fields.bzfield = fieldgroup.member[i][j].bfield.member[2];
        }
    }
}
/// 演化磁场，并damping
void __Simbox::advancest1(const int& rank)
{
    if(disable_em)
    {
        return;
    }
    //cout<<"initing flag, em will be executed "<<endl;
    init_flag();
    fieldgroup.advancepre(rank, delta_t);
}

/// 演化电场damping
void __Simbox::advancest2(const int& rank)
{
    if(disable_em)
    {
        return;
    }
    fieldgroup.update_efield(delta_t);
}

/// 演化磁场反damping，拷贝场到cell，场weight到ptcls
void __Simbox::advancest3(const int& rank)
{
    if(disable_em)
    {
        return;
    }
#ifdef __Debug
    cout<<"field advancenex "<<endl;
#endif
    fieldgroup.advancenex(rank, delta_t);
#ifdef __Debug
    cout<<" field copy to cell "<<endl;
#endif
    field_copy_to_cell();
#ifdef __Debug
    cout<<" weighting to ptcls "<<endl;
#endif
    weight_to_ptcls();
#ifdef __Debug
    cout<<" weighting over "<<endl;
#endif
}

/// 推进粒子
void __Simbox::advancest4(const int& rank, const int& n)
{
    //cout<<"size "<<cellgroup.member[30][30].plist[0].size()<<endl;
    // 这两个不能交换，否则粒子被压栈到ｂｄ中就不能辐射了
    
    for(int i = 0; i < 4; i ++)
    {
        if(monte_flag[i]!=0)
        {
#ifdef __Debug
            cout<<rank<<" begin monte_carlo "<<i<<endl;
#endif
            if(i == 0 || i == 2)
            {
                monte_carlo(cellgroup, i, delta_t, tables1d, tables2d);
            }
            else{
                monte_carlo(cellgroup, i, delta_t, bh_tables1d, bh_tables2d);
            }
#ifdef __Debug
            cout<<rank<<" over monte_carlo "<<i<<endl;
#endif
        }
    }
    
    cellgroup.advance(delta_t, crossing_bd_ptcls, rank, n);
}
// 交换数据后只后要进行migration

/// 进行migration，将粒子weight到grid上，拷贝到field中
void __Simbox::advancest5(const int& rank)
{
    migration(rank);
    /// 必须开启,因为密度更新也在此过程中.
    cellgroup.weight_to_grid(rank);
    //cellgroup.weight_to_gridc();
    //cellgroup.weight_to_gridb();
    //cout<<rank<<" w t g over "<<endl;
}


__Vect2<int> __Simbox::coords_to_cell(const __Vect2<double>& xycoords)
{
    __Vect2<int> cell_position(-1, -1);
    __Vect2<double> temp = xycoords - xymin;//cellgroup.member[0][0].xymin;
    if(temp.member[0] >= 0 && temp.member[1] >= 0)
    {
        cell_position.member[0] = temp.member[0] / dxy.member[0];
        cell_position.member[1] = temp.member[1] / dxy.member[1];
        if(cell_position.member[0] < cellgroup.member.size() &&\
                cell_position.member[0] >= 0 && \
                cell_position.member[1] < cellgroup.member[0].size() &&\
                cell_position.member[1] >= 0)
        {
            return cell_position;
        }
        else{
            cell_position = __Vect2<int> (-1, -1);
        }
    }
    else{
        cell_position = __Vect2<int> (-1, -1);
    }
    return cell_position;
}


void __Simbox::migration(const int& myrank)
{

    for(int i = 0; i < cellgroup.member[0].size(); i ++)
    {
        for(int j = 0; j < species; j ++)
        {
            cellgroup.member[0][i].plist[j].clear();
            cellgroup.member[cellgroup.member.size() - 1][i].plist[j].clear();
        }
    }
    for(int i = 0; i < cellgroup.member.size(); i ++)
    {
        for(int j = 0; j < species; j ++)
        {
            cellgroup.member[i][0].plist[j].clear();
            cellgroup.member[i][cellgroup.member[i].size() - 1].plist[j].clear();
        }
    }

    for(int i = 0; i < this->species; i++)
    {
        int total_size = crossing_bd_ptcls[i].size();
        for(int j = total_size - 1; j >= 0; j --)
        {
            __Ptcls& thisptcl = crossing_bd_ptcls[i][j];
            //if((xymax >= thisptcl.position) && (thisptcl.position >= xymin))
            //if ((thisptcl.position.member[0] <= this->xymax.member[0]) && (thisptcl.position.member[0] >= \
                this->xymin.member[0]) && (thisptcl.position.member[1] <= xymax.member[1]) && (thisptcl.position.\
                    member[1] >= xymin.member[0]))
            //{
                __Vect2<int> cell_position = coords_to_cell(thisptcl.position);
                if(cell_position.member[0] >= 0 && cell_position.member[1] >= 0)
                {
                    this->cellgroup.member[cell_position.member[0]]\
                    [cell_position.member[1]].push(thisptcl);
                }
            //}
            else
            {
                /// particles in this domain will be deleted
                __Ptcls::total_num --;
            }
        }
        clear_vec(crossing_bd_ptcls[i]);
    }
}

int __Simbox::read_tables2d()
{
    //cout<<"reading 2d size"<<filenames2d.size()<<endl;
    vector<vector<double> > tabletemp;
    for(int i = 0; i < filenames2d.size(); i ++)
    {
        //cout<<"file name size "<<filenames2d.size()<<endl;
        cout<<"reading "<<filenames2d[i]<<endl;
        ifstream in(filenames2d[i], ios::in);
        double temp;
        string line;
        vector<double> row;
        while(getline(in, line)){
            istringstream istr(line);
            while(istr>>temp){
                //cout<<"2d reading "<<temp<<endl;
                row.push_back(temp);
            }
            tabletemp.push_back(row);
            istr.clear();
            line.clear();
            row.clear();
        }
        if(i > 2)
        {
            bh_tables2d.push_back(tabletemp);
        }
        else{
            tables2d.push_back(tabletemp);
        }
        //cout<<"temp size "<<tabletemp.size()<<endl;
        clear_vec(tabletemp);
        //tabletemp.clear();
    }
    return 0;
}



int __Simbox::read_tables1d()
{
    double temp;
    for(int i = 0; i < filenames1d.size(); i ++)
    {
        cout<<"reading "<<filenames1d[i]<<endl;
        ifstream in;
        in.open(filenames1d[i], ios::in);
        vector<double> tb;
        while(in>>temp)  
        {
            tb.push_back(temp);
        }
        if(i > 6)
        {
            bh_tables1d.push_back(tb);
        }
        else{
            tables1d.push_back(tb);
        }
        clear_vec(tb);
    }
    return 0;
}



void __Simbox::init()
{
    vector<__Ptcls> tempbundle;
    for(int i = 0; i < species; i ++)
    {
        crossing_bd_ptcls.push_back(tempbundle);
    }
    //cout<<"initing "<<crossing_bd_ptcls[0].size()<<endl;
    //cout<<"creating cell group"<<endl;
    __Cellgroup mygroup(xygridmin, xygridmax, xymin, xymax, xycell, species, ppc, spptcl, \
            delta_t, position, maxid);
    mygroup.dxy = dxy;
    cellgroup = mygroup;
    cellgroup.disable_em = disable_em;
    //cout<<"cell group initing"<<endl;
    cellgroup.init();
    cellgroup.set_bd(bd_type);
    //cout<<"initing field"<<endl;
    __Fieldgroup myfield(xycell.member[0] + 1, xycell.member[1] + 1, \
            dxy, bd_type, startindex, damplen);
    fieldgroup = myfield;
    fieldgroup.init();
 
    /***
    cellgroup.dxy = dxy;
    cellgroup.disable_em = disable_em;
    cout<<"cell group initing"<<endl;
    cellgroup.init();
    cout<<"cellgroup init over "<<endl;
    cellgroup.set_bd(bd_type);
    //cout<<"initing field"<<endl;
    fieldgroup.init();
    ***/


    if(monte_flag[0] + monte_flag[1] + monte_flag[2] + monte_flag[3] != 0)
    {
        if( read_tables1d() != 0)
        {
            cout<<"read tables1d failed, monte_carlo will not work!"<<endl;
            for(int i = 0; i < 4; i ++)
            {
                monte_flag[i] = 0;
            }
        }
        if(read_tables2d() != 0)
        {
            cout<<"read tables2d failed, monte_carlo will not work!"<<endl;
            for(int i = 0; i < 4; i ++)
            {
                monte_flag[i] = 0;
            }
        }
    }
}


ostream& operator<<(ostream &out, const __Simbox& mybox){
    out<<" simname: "<<mybox.simname<<endl<<" date: "<<mybox.simdate<<endl\
    <<" simx: "<<mybox.xrange<<endl<<" simy: "<<mybox.yrange<<endl<<" xycell: "<<\
    mybox.xycell<<endl<<" xygridmin: "<<mybox.xygridmin<<endl<<" xygridmax: "<<mybox.xygridmax\
    <<endl<<" time: "<<mybox.time;
	return out;
}


// direction 0: left & right, 1: down & up
void __Simbox::copy_to_pline(__Pline& plinel, __Pline& pliner, int direction)
{
    __Pdata pdata;
    for(int i = 0; i < crossing_bd_ptcls.size(); i ++)
    {
        for(int j = 0; j < crossing_bd_ptcls[i].size(); j ++)
        {
            __Ptcls& ptcls = crossing_bd_ptcls[i][j];
            if(direction == 0)
            {
                //if(ptcls.position.member[0] < xymin.member[0] )
                if(ptcls.position.member[0] < xymin.member[0] + dxy.member[0])
                {
                    copy_to_pdata(ptcls, pdata);
                    plinel.push_back(pdata);
                }
                else if(ptcls.position.member[0] > xymax.member[0] - dxy.member[0])
                {
                    copy_to_pdata(ptcls, pdata);
                    pliner.push_back(pdata);
                }
                
                // 这里由于越界例子会在migration时删除因此，我不需要把它们从crossing_bd中删除
            }
            if(direction == 1)
            {
                if(ptcls.position.member[1] < xymin.member[1] + dxy.member[1])
                {
                    copy_to_pdata(ptcls, pdata);
                    plinel.push_back(pdata);
                }
                else if(ptcls.position.member[1] > xymax.member[1] - dxy.member[1])
                {
                    copy_to_pdata(ptcls, pdata);
                    pliner.push_back(pdata);
                }
                // 这里由于越界例子会在migration时删除因此，我不需要把它们从crossing_bd中删除
            }
        }
    }
    if(direction == 0)
    {
        for(int i = 0; i < cellgroup.member[0].size(); i ++)
        {
            for(int j = 0; j < species; j ++)
            {
                for(int k = 0; k < cellgroup.member[1][i].plist[j].size(); k ++)
                {
                    copy_to_pdata(cellgroup.member[1][i].plist[j].member[k], pdata);
                    plinel.push_back(pdata);
                }
                for(int k = 0; k < cellgroup.member[cellgroup.member.size() - 2][i].plist[j].size(); k ++)
                {
                    copy_to_pdata(cellgroup.member[cellgroup.member.size() - 2][i].plist[j].member[k], pdata);
                    pliner.push_back(pdata);
                }
            }
        }
    }
    if(direction == 1)
    {
        for(int i = 0; i < cellgroup.member.size(); i ++)
        {
            for(int j = 0; j < species; j ++)
            {
                for(int k = 0; k < cellgroup.member[i][1].plist[j].size(); k ++)
                {
                    copy_to_pdata(cellgroup.member[i][1].plist[j].member[k], pdata);
                    plinel.push_back(pdata);
                }
                for(int k = 0; k < cellgroup.member[i][cellgroup.member.size() - 2].plist[j].size(); k ++)
                {
                    copy_to_pdata(cellgroup.member[i][cellgroup.member.size() - 2].plist[j].member[k], pdata);
                    pliner.push_back(pdata);
                }
            }
        }
    }

}



void __Simbox::pline_copy_to(__Pline& plinel, __Pline& pliner, int direction)
{
    /***
    __Spinfo spinfo[4] = {{0, e_mass, - e_charge}, {1, 1837.0 * e_mass * ATOM_NUM, \
        e_charge * ION_CHARGE}, {2, e_mass, e_charge}, {3, 0, 0}};
    ***/
    // 这个其实跟方向就没有关系了,因为跟网格暂时没有关系
    if(plinel.size() + pliner.size() == 0)
    {
        return;
    }
    __Ptcls ptcls;
    if(direction == 0)
    {
        for(int i = 0; i < cellgroup.member[0].size(); i ++)
        {
            for(int j = 0; j < species; j ++)
            {
                cellgroup.member[0][i].plist[j].clear();
                cellgroup.member[cellgroup.member.size() - 1][i].plist[j].clear();
            }
        }
        for(int i = 0; i < plinel.size(); i ++)
        {
            copy_to_ptcl(plinel[i], ptcls);
            crossing_bd_ptcls[ptcls.species].push_back(ptcls);
            ptcls.total_id ++;
        }
        for(int i = 0; i < pliner.size(); i ++)
        {
            copy_to_ptcl(pliner[i], ptcls);
            crossing_bd_ptcls[ptcls.species].push_back(ptcls);
            ptcls.total_id ++;
        }
    }
    
    
    if(direction == 1)
    {
        for(int i = 0; i < cellgroup.member.size(); i ++)
        {
            for(int j = 0; j < species; j ++)
            {
                cellgroup.member[i][0].plist[j].clear();
                cellgroup.member[i][cellgroup.member[0].size() - 1].plist[j].clear();
            }
        }
        for(int i = 0; i < plinel.size(); i ++)
        {
            copy_to_ptcl(plinel[i], ptcls);
            crossing_bd_ptcls[ptcls.species].push_back(ptcls);
            ptcls.total_id ++;
        }
        for(int i = 0; i < pliner.size(); i ++)
        {
            copy_to_ptcl(pliner[i], ptcls);
            crossing_bd_ptcls[ptcls.species].push_back(ptcls);
            ptcls.total_id ++;
        }
    }
}


void __Simbox::fline_copy_to(__Fline& flinel, __Fline& fliner, int direction)
{
    //接收两列，总是先左后右，线下后上。
    //cout<<"fline copy to"<<endl;
    if(direction == 0)
    {
        for(int i = 0; i < flinel.size() / 2; i ++)
        {
            fieldgroup.member[0][i] = flinel[i];
            /***
            __Vect3<double> myjfield = fieldgroup.member[1][i].jfield;
            double rho = fieldgroup.member[1][i].rho;
            ***/
            fieldgroup.member[1][i] = flinel[i + flinel.size() / 2];
            /***
            fieldgroup.member[1][i].jfield = myjfield + fieldgroup.member[1][i].jfield;
            fieldgroup.member[1][i].rho += rho;
            ***/
        }
        int ind = fieldgroup.member.size();
        for(int i = 0; i < fliner.size() / 2; i ++)
        {
            fieldgroup.member[ind - 1][i] = fliner[i + fliner.size() / 2];
            /***
            __Vect3<double> myjfield = fieldgroup.member[ind - 2][i].jfield;
            double rho = fieldgroup.member[ind - 2][i].rho;
            ***/
            fieldgroup.member[ind - 2][i] = fliner[i];
            /***
            fieldgroup.member[ind - 2][i].jfield = myjfield + fieldgroup.member[ind - 2][i].jfield;
            fieldgroup.member[ind - 2][i].rho += rho;
            ***/
        }
    }
    if(direction == 1)
    {
        for(int i = 0; i < flinel.size() / 2; i ++)
        {
            fieldgroup.member[i][0] = flinel[i];
            /***
            __Vect3<double> myjfield = fieldgroup.member[i][1].jfield;
            double rho = fieldgroup.member[i][1].rho;
            ***/
            fieldgroup.member[i][1] = flinel[i + flinel.size() / 2];
            /***
            fieldgroup.member[i][1].jfield = myjfield + fieldgroup.member[i][1].jfield;
            fieldgroup.member[i][1].rho += rho;
            ***/
            
        }
        int ind = fieldgroup.member[0].size();
        for(int i = 0; i < fliner.size() / 2; i ++)
        {
            fieldgroup.member[i][ind - 1] = fliner[i + fliner.size() / 2];
            /***
            __Vect3<double> myjfield = fieldgroup.member[i][ind - 2].jfield;
            double rho = fieldgroup.member[i][ind - 2].rho;
            ***/
            fieldgroup.member[i][ind - 2] = fliner[i];
            /***
            fieldgroup.member[i][ind - 2].jfield = myjfield + fieldgroup.member[i][ind - 2].jfield;
            fieldgroup.member[i][ind - 2].rho += rho;
            ***/

        }
    }
    //cout<<"fline copy to over "<<endl;
}
void __Simbox::copy_to_fline(__Fline& flinel, __Fline& fliner, int direction)
{
    //cout<<"copy to fline"<<endl;
    __Fdata fdatal, fdatar;
    if(direction == 0)
    {
        for(int i = 0; i < fieldgroup.member[1].size(); i ++)
        {
            copy_to_fdata(fieldgroup.member[1][i], fdatal);
            flinel.push_back(fdatal);
            copy_to_fdata(fieldgroup.member[fieldgroup.member.size() - 3][i], fdatar);
            fliner.push_back(fdatar);
        }
        for(int i = 0; i < fieldgroup.member[1].size(); i ++)
        {
            copy_to_fdata(fieldgroup.member[2][i], fdatal);
            flinel.push_back(fdatal);
            copy_to_fdata(fieldgroup.member[fieldgroup.member.size() - 2][i], fdatar);
            fliner.push_back(fdatar);
        }
    }
    if(direction == 1)
    {
        for(int i = 0; i < fieldgroup.member.size(); i ++)
        {
            copy_to_fdata(fieldgroup.member[i][1], fdatal);
            flinel.push_back(fdatal);
            copy_to_fdata(fieldgroup.member[i][fieldgroup.member[0].size() - 3], fdatar);
            fliner.push_back(fdatar);
        }
        for(int i = 0; i < fieldgroup.member.size(); i ++)
        {
            copy_to_fdata(fieldgroup.member[i][2], fdatal);
            flinel.push_back(fdatal);
            copy_to_fdata(fieldgroup.member[i][fieldgroup.member[0].size() - 2], fdatar);
            fliner.push_back(fdatar);
        }
    }
    //cout<<"copy to fline over"<<endl;
    
}

void __Simbox::copy_to_jline(__Jline& jlinel, __Jline& jliner, int direction)
{
    //cout<<"start copy to jline "<<endl;
    __Jdata jdatal, jdatar;
    int xsize = cellgroup.jfield.size(), ysize = cellgroup.jfield[0].size();
    if(direction == 0)
    {
        for(int i = 0; i < cellgroup.jfield[0].size(); i ++)
        {
            //cout<<" copying "<<i<<" xsize - 3 "<<xsize - 3<<endl;
            copy_to_jdata(cellgroup.jfield[1][i], cellgroup.rhofield[1][i], jdatal);
            jlinel.push_back(jdatal);
            copy_to_jdata(cellgroup.jfield[xsize - 3][i], cellgroup.rhofield[xsize - 3][i], jdatar);
            jliner.push_back(jdatar);
        }
        //cout<<"copy over"<<endl;
        for(int i = 0; i < cellgroup.jfield[0].size(); i ++)
        {
            //cout<<" copying right "<<i<<" xsize - 3 "<<xsize - 3<<endl;
            copy_to_jdata(cellgroup.jfield[2][i], cellgroup.rhofield[2][i], jdatal);
            jlinel.push_back(jdatal);
            copy_to_jdata(cellgroup.jfield[xsize - 2][i], cellgroup.rhofield[xsize - 2][i], jdatar);
            jliner.push_back(jdatar);
        }
    }
    if(direction == 1)
    {
        for(int i = 0; i < cellgroup.jfield.size(); i ++)
        {
            /// 左和下相同，右和上相同
            copy_to_jdata(cellgroup.jfield[i][1], cellgroup.rhofield[i][1], jdatal);
            jlinel.push_back(jdatal);
            copy_to_jdata(cellgroup.jfield[i][ysize - 3], cellgroup.rhofield[i][ysize - 3], jdatar);
            jliner.push_back(jdatar);
        }
        for(int i = 0; i < cellgroup.jfield.size(); i ++)
        {
            /// 左和下相同，右和上相同
            copy_to_jdata(cellgroup.jfield[i][2], cellgroup.rhofield[i][2], jdatal);
            jlinel.push_back(jdatal);
            copy_to_jdata(cellgroup.jfield[i][ysize - 2], cellgroup.rhofield[i][ysize - 2], jdatar);
            jliner.push_back(jdatar);
        }
    }
    //cout<<" jline over "<<endl;

}


/// 将信息储存在边缘，最后叠加到１上
void __Simbox::jline_copy_to(__Jline& jlinel, __Jline& jliner, int direction)
{
    //cout<<"start jline copy to "<<endl;
    int posl = jlinel.size() / 2, posr = jliner.size() / 2;
    int xsize = cellgroup.jfield.size(), ysize = cellgroup.jfield[0].size();
    if(direction == 0)
    {
        for(int i = 0; i < posl; i ++)
        {
            //cout<<"ladding to "<<i<<endl;
            for(int j = 0; j < 3; j ++)
            {
                cellgroup.jfield[0][i].member[j] = jlinel[i].data[j];
                cellgroup.jfield[1][i].member[j] += jlinel[i + posl].data[j];
            }
            cellgroup.rhofield[0][i] = jlinel[i].data[3];
            cellgroup.rhofield[1][i] += jlinel[i + posl].data[3];
        }
        for(int i = 0; i < posr; i ++)
        {
            //cout<<"radding to "<<i<<endl;
            for(int j = 0; j < 3; j ++)
            {
                cellgroup.jfield[xsize - 2][i].member[j] += jliner[i].data[j];
                cellgroup.jfield[xsize - 1][i].member[j] = jliner[i + posr].data[j];
            }
            cellgroup.rhofield[xsize - 2][i] += jliner[i].data[3];
            cellgroup.rhofield[xsize - 1][i] = jliner[i + posr].data[3];
        }
    }
    if(direction == 1)
    {
        for(int i = 0; i < jlinel.size() / 2; i ++)
        {
            for(int j = 0; j < 3; j ++)
            {
                cellgroup.jfield[i][0].member[j] = jlinel[i].data[j];
                cellgroup.jfield[i][1].member[j] += jlinel[i + jlinel.size() / 2].data[j];
            }
            cellgroup.rhofield[i][0] = jlinel[i].data[3];
            cellgroup.rhofield[i][1] += jlinel[i + jlinel.size() / 2].data[3];
        }
        for(int i = 0; i < jliner.size() / 2; i ++)
        {
            for(int j = 0; j < 3; j ++)
            {
                cellgroup.jfield[i][ysize - 2].member[j] += jliner[i].data[j];
                cellgroup.jfield[i][ysize - 1].member[j] = jliner[i + jliner.size() / 2].data[j];
            }
            cellgroup.rhofield[i][ysize - 2] += jliner[i].data[3];
            cellgroup.rhofield[i][ysize - 1] = jliner[i + jliner.size() / 2].data[3];
        }
    }
}






/***

/// 采用T.D.Pointon/CPC 179(2008) 535-544 算法
void __Simbox::weight_to_ptcls()
{
    __Vect3<double> zeros(0.0, 0.0, 0.0);

    for(int i = 1; i < cellgroup.member.size() - 1; i ++)
    {
        for(int j = 1; j < cellgroup.member[1].size() - 1; j ++)
        {
            __Vect2<double>& cellmin = cellgroup.member[i][j].xymin;
            for(int k = 0; k < cellgroup.member[i][j].plist.size(); k ++)
            {
                for(int l = 0; l < cellgroup.member[i][j].plist[k].size(); l ++)
                {
                    __Ptcls& myptcl = cellgroup.member[i][j].plist[k].member[l];
                    myptcl.efield = zeros;
                    myptcl.bfield = zeros;
                    /// 对相邻的三个cell中的场进行插值到粒子
                    /// 这里采用相对坐标local
                    int indx, indy;
                    __Vect2<double> ptcls_ratio(myptcl.position - cellmin);
                    ptcls_ratio.member[0] /= dxy.member[0];
                    ptcls_ratio.member[1] /= dxy.member[1];
                    indx = i, indy = j;
                    if(ptcls_ratio.member[0] > 0.5)
                    {
                        indx = i + 1;
                        ptcls_ratio.member[0] = 1.0 - ptcls_ratio.member[0];
                    }
                    if(ptcls_ratio.member[1] > 0.5)
                    {
                        indy = j + 1;
                        ptcls_ratio.member[1] = 1.0 - ptcls_ratio.member[1];
                    }
                    double wx[3], wy[3], ux[2], uy[2];
                    ux[0] = ptcls_ratio.member[0], ux[1] = 1.0 - ux[0];
                    uy[0] = ptcls_ratio.member[1], uy[1] = 1.0 - uy[0];
                    wx[0] = 0.5 * ux[0] * ux[0], wx[2] = 0.5 * ux[1] * ux[1], wx[1] = 1.0 - wx[0] - wx[2];
                    wy[0] = 0.5 * uy[0] * uy[0], wy[2] = 0.5 * uy[1] * uy[1], wy[1] = 1.0 - wy[0] - wy[2];
                    
                    /// 3 grid near the current cell, since field are all locate on the center
                    for(int yj = -1; yj < 2; yj ++)
                    {
                        myptcl.efield.member[0] += (fieldgroup.member[indx - 1][indy + yj].efield.member[0] * \
                                ux[0] + fieldgroup.member[indx][indy + yj].efield.member[0] * ux[1]) * wy[yj + 1];
                        myptcl.bfield.member[0] += (fieldgroup.member[indx - 1][indy+yj].bfield.member[0] * \
                                ux[0] + fieldgroup.member[indx][indy + yj].bfield.member[0] * ux[1]) * wy[yj + 1];
                    }
                    
                    for(int xi = -1; xi < 2; xi ++)
                    {
                        myptcl.efield.member[1] += (fieldgroup.member[indx+xi][indy - 1].efield.member[1] * \
                                uy[0] + fieldgroup.member[indx + xi][indy].efield.member[1] * uy[1]) * wx[xi + 1];
                        myptcl.bfield.member[1] += (fieldgroup.member[indx+xi][indy - 1].bfield.member[1] * \
                                uy[0] + fieldgroup.member[indx + xi][indy].bfield.member[1] * uy[1]) * wx[xi + 1];
                    }
                    
                    for(int xi = -1; xi < 2; xi ++)
                    {
                        for(int yj = -1; yj < 2; yj ++)
                        {
                            myptcl.efield.member[2] += fieldgroup.member[indx+xi][indy+yj].efield.member[2] * \
                                                    wx[xi + 1] * wy[yj + 1];
                            myptcl.bfield.member[2] += fieldgroup.member[indx+xi][indy+yj].bfield.member[2] * \
                                                    wx[xi + 1] * wy[yj + 1];
                        }
                    }
                }
            }
        }
    }
}
***/
/// 采用如同EPOCH中的粒子三角形，场均匀算法
void __Simbox::weight_to_ptcls()
{
    __Vect3<double> zeros(0.0, 0.0, 0.0);
    __Vect2<double> xhalf(dxy.member[0] * 0.5, 0);
    __Vect2<double> yhalf = dxy * 0.5 - xhalf;
    double frac = _frac2;
    for(int i = 1; i < cellgroup.member.size() - 1; i ++)
    {
        for(int j = 1; j < cellgroup.member[1].size() - 1; j ++)
        {
            int indx = i, indy = j;
            __Vect2<double>& cellmin = cellgroup.member[i][j].xymin;
            for(int k = 0; k < cellgroup.member[i][j].plist.size(); k ++)
            {
                for(int l = 0; l < cellgroup.member[i][j].plist[k].size(); l ++)
                {
                    __Ptcls& myptcl = cellgroup.member[i][j].plist[k].member[l];
                    if(std::isnan(myptcl.energy)) cout<<myptcl<<" nan energy"<<endl;
                    myptcl.efield = zeros;
                    myptcl.bfield = zeros;
                    /// 对相邻的三个cell中的场进行插值到粒子
                    /// 这里采用相对坐标local
                    
                    ///已经验证，没有错误
                    indx = floor((myptcl.position.member[0] + dxy.member[0] * 0.5 - \
                                cellmin.member[0]) / dxy.member[0]) + i;
                    indy = floor((myptcl.position.member[1] + dxy.member[1] * 0.5 - \
                                cellmin.member[1]) / dxy.member[1]) + j;
                    __Vect2<double> ptcls_ratio((cellgroup.member[indx][indy].xymin - myptcl.position) / dxy);
                    __Vect2<double> ptcls_ratio_center;
                    /// 对中间点进行shift
                    ptcls_ratio_center.member[0] = (ptcls_ratio.member[0] > 0)?(0.5 - ptcls_ratio.member[0]):\
                                                   (0.5 + ptcls_ratio.member[0]);
                    ptcls_ratio_center.member[1] = (ptcls_ratio.member[1] > 0)?(0.5 - ptcls_ratio.member[1]):\
                                                   (0.5 + ptcls_ratio.member[1]);
                    double tempx0[5], tempy0[5], tempx1[5], tempy1[5];
                    //ptcls_ratio_center = ptcls_ratio;
                    grid_wxy(ptcls_ratio.member[0], ptcls_ratio.member[1], tempx0, tempy0);
                    grid_wxy(ptcls_ratio_center.member[0], ptcls_ratio_center.member[1], tempx1, tempy1);



                    /// 3 grid near the current cell, since field are all locate on the center
                    for(int xi = _xmin; xi <= _xmax; xi ++)
                    {
                        if(indx + xi < 1 || indx + xi > fieldgroup.member.size() - 1)
                        {
                            continue;
                        }
                        for(int yj = _ymin; yj <= _ymax; yj ++)
                        {
                            if(indy + yj < 1 || indy + yj > fieldgroup.member[0].size() - 1)
                            {
                                continue;
                            }
                            
                            /// use square area weighting
                            /***
                            __Vect2<double> local = xhalf * (2.0 * xi) + yhalf * (2.0 * yj);
                            myptcl.efield.member[0] += fieldgroup.member[indx+xi][indy+yj].efield.member[0] \
                                                       * weightfunc(myptcl.position, cellmin + local + xhalf, dxy);
                            myptcl.efield.member[1] += fieldgroup.member[indx+xi][indy+yj].efield.member[1] \
                                                       * weightfunc(myptcl.position, cellmin + local + yhalf, dxy);
                            myptcl.efield.member[2] += fieldgroup.member[indx+xi][indy+yj].efield.member[2] \
                                                       * weightfunc(myptcl.position, cellmin + local + xhalf + yhalf, dxy);
                            myptcl.bfield.member[0] += fieldgroup.member[indx+xi][indy+yj].bfield.member[0] \
                                                       * weightfunc(myptcl.position, cellmin + local + xhalf, dxy);
                            myptcl.bfield.member[1] += fieldgroup.member[indx+xi][indy+yj].bfield.member[1] \
                                                       * weightfunc(myptcl.position, cellmin + local + yhalf, dxy);
                            myptcl.bfield.member[2] += fieldgroup.member[indx+xi][indy+yj].bfield.member[2] \
                                                       * weightfunc(myptcl.position, cellmin + local + xhalf + yhalf, dxy);
                            ***/

                            myptcl.efield.member[0] += fieldgroup.member[indx+xi][indy+yj].efield.member[0] \
                                                        * tempx1[xi + 2] * tempy0[yj + 2] * frac;
                            myptcl.efield.member[1] += fieldgroup.member[indx+xi][indy+yj].efield.member[1] \
                                                        * tempx1[xi + 2] * tempy0[yj + 2] * frac;
                            myptcl.efield.member[2] += fieldgroup.member[indx+xi][indy+yj].efield.member[2] \
                                                        * tempx1[xi + 2] * tempy1[yj + 2] * frac;
                            myptcl.bfield.member[0] += fieldgroup.member[indx+xi][indy+yj].bfield.member[0] \
                                                        * tempx1[xi + 2] * tempy0[yj + 2] * frac;
                            myptcl.bfield.member[1] += fieldgroup.member[indx+xi][indy+yj].bfield.member[1] \
                                                        * tempx0[xi + 2] * tempy1[yj + 2] * frac;
                            myptcl.bfield.member[2] += fieldgroup.member[indx+xi][indy+yj].bfield.member[2] \
                                                        * tempx1[xi + 2] * tempy1[yj + 2] * frac;
                        }
                    }
                }
            }
        }
    }
}



void __Simbox::wide_merge(const int& sp, const int& cellnum)
{
    for(int cidx = cellnum; cidx < cellgroup.member.size() - cellnum; cidx ++)
    {
        for(int cidy = cellnum; cidy < cellgroup.member[0].size() - cellnum; cidy++)
        {
            __Ptclslist& current = cellgroup.member[cidx][cidy].plist[sp];
            int nid = 0;
            __Vect3<double> dvv;
            while(nid < current.member.size())
            {
                for(int xid = 0; xid < cellnum; xid ++)
                {
                    for(int yid = 0; yid < cellnum; yid ++)
                    {
                        __Ptclslist& mlist = cellgroup.member[xid][yid].plist[sp];
                        for(auto p = mlist.member.begin(); p != mlist.member.end(); p ++)
                        {
                            dvv = current.member[nid].velocity - p->velocity;
                            double delta_energy = abs(current.member[nid].energy - p->energy);
                            double delta_velocity = Mod(dvv);
                            /// 能量上加上1ev,防止粒子的能量为0, 速度上加1m/s,防止粒子完全静止
                            if(delta_energy / (current.member[nid].energy + 1e-19) <= delta_e && \
                                    (delta_velocity / (Mod(current.member[nid].velocity) + 1.0) < delta_v))
                            //if(delta_energy / member[nid].energy <= delta_e && delta_energy / p->energy <= delta_e && \
                                    (delta_velocity / Mod(member[nid].velocity) < delta_v && Mod(dvv) / \
                                        Mod(p->velocity) < delta_v))
                            {
                                current.member[nid].energy = (current.member[nid].energy * \
                                        current.member[nid].weight + p->energy * p->weight);
                                current.member[nid].velocity = (current.member[nid].velocity * \
                                        current.member[nid].weight + p->velocity * p->weight);
                                current.member[nid].position = (current.member[nid].position + p->position) * 0.5;
                                current.member[nid].weight += p->weight;
                                current.member[nid].energy /= current.member[nid].weight;
                                current.member[nid].velocity = current.member[nid].velocity / current.member[nid].weight;
                                /// 删除被merge的粒子 
                                p = mlist.member.erase(p);
                                p --;
                            }
                        }
                        mlist.squeeze();
                        nid ++;
                    }
                    //cout<<" "<<member.size()<<endl;
                    /// 进行压缩空间
                }
            }
        }
    }
}
                       




#endif
